#### Code to produce results reported in the main text  ####
#### Table 1, Figure 3, Table 2, and Figure 2 (in that order) ####
## Author: Jeong Hyun Kim
## Last updated: 12/28/2021

library(foreign)
library(tidyverse)
library(broom)
library(stargazer)

#### Study 1: Loading Data and Coding variables ####
csid <- read.csv("survey_combined.csv", stringsAsFactors = FALSE)
csid$gov_support <- ifelse(csid$BQ14=="1", 1, 0)
# Coding young1: those in 20s and 30s
csid$young1 <- ifelse(csid$BQ3 < 40, 1, 0)
# Alternative operationalization of young variable
csid$young2 <- ifelse(csid$BQ3 < 40 & csid$BQ3 > 25, 1, 0)
# Alternative operationalization of young variable
csid$young3 <- ifelse(csid$BQ3 < 39 & csid$BQ3 > 22, 1, 0)  
names(csid)[3] <- "region"
csid$region <- as.factor(csid$region)
names(csid)[4] <- "age"
names(csid)[5] <- "education"
names(csid)[6] <- "occupation"
csid$occupation <- as.factor(csid$occupation)
names(csid)[7] <- "job_status"
csid[which(is.na(csid$job_status)), "job_status"] <- "other"
csid[which(csid$job_status=="#NULL!"),"job_status"] <- "other" # unemployed; retired; student
csid[which(csid$job_status=="1"),"job_status"] <- "regular"
csid[which(csid$job_status=="2"),"job_status"] <- "temp"
csid$temp_worker <- ifelse(csid$job_status=="temp", 1, 0)


names(csid)[8] <- "marital_status"
names(csid)[9] <- "children"
names(csid)[10] <- "religion"
names(csid)[11] <- "income"
names(csid)[12] <- "ideology"
names(csid)[13] <- "partyID"

# For BQ13, 99 is coded as don't know.
csid[which(csid$ideology=="99"),"ideology"] <- NA

#### Table 1, Panel (a) ####
mod1.1 <- lm(gender_quota_att ~ treated +  as.factor(wave), csid)
summary(mod1.1)
nrow(mod1.1$model)

mod1.2 <- lm(gender_quota_att ~ treated +  as.factor(wave) + young1 
             + region + education + income + as.factor(occupation) 
             + ideology + as.factor(partyID) 
             + econ_bad + as.factor(religion) +  as.factor(marital_status) 
             + as.factor(children), csid)
summary(mod1.2)
nrow(mod1.2$model)

mod1.3 <- lm(gender_quota_att ~ treated*young1 +  as.factor(wave) , csid)
summary(mod1.3)
nrow(mod1.3$model)

mod1.4 <- lm(gender_quota_att ~ treated*young1 +  as.factor(wave) +   region
             + education + income + as.factor(occupation) + ideology + as.factor(partyID)
             + econ_bad + as.factor(religion) + as.factor(marital_status) 
             + as.factor(children), csid)

summary(mod1.4)
nrow(mod1.4$model)


#### Study 2: Loading Data and Coding Variables ####
df <- read.spss("survey2021.sav", to.data.frame = TRUE, use.value.labels = FALSE)

# Create a variable that denotes respondents who passed manipulation check:
df <- df %>% mutate(T1_pass = case_when(K0 == 1 & K1==1  ~ 1,
                                        K0 == 1 & K1==2 ~ 0),
                    T2_pass = case_when(K0 == 2 & K2 ==1 ~ 1,
                                        K0 == 2 & K2 == 2 ~ 0))
table(df[which(df$K0==1), "T1_pass"])
table(df[which(df$K0==2), "T2_pass"])

# Reduce the sample to those who passed the manipulation check:
df <- df[which(df$T1_pass==1 | df$T2_pass==1), ]

# code the respondents who received "T1" as treated:
df <- df %>% mutate(treated = case_when(K0 == 1 ~ 1,
                                        K0 == 2 ~ 0))

df$young1 <- ifelse(df$BQ3 < 40, 1, 0)
# Alternative operationalization of young variable
df$young2 <- ifelse(df$BQ3 < 40 & df$BQ3 > 25, 1, 0)
# Alternative operationalization of young variable
df$young3 <- ifelse(df$BQ3 < 41 & df$BQ3 > 24, 1, 0)  
df$young4 <- ifelse(df$BQ3 < 35 & df$BQ3 > 19, 1, 0)  # 20-34

df<-  df %>% rename(age = BQ3,
                    education = BQ4,
                    occupation = BQ5,
                    job_status = BQ5_3,
                    region = BQ1,
                    income = BQ12_1,
                    ideology = BQ13,
                    religion = BQ8,
                    marital_status = BQ6,
                    gender_quota_att = K3_1,
                    equal_pay = K3_2,
                    corporate_quota = K3_3
)
df <- df %>% rename(male_status = K6,
                    female_status = K7)

df <- df %>% mutate(children = case_when(ABQ7 == 1 ~ 1,
                                         ABQ7 == 2 ~ 0),
                    econ_bad = ifelse(K9==1| K9==2, 1, 0))
## gender norms
# K5_1: 가사노동 (lower value indicates more related to men)
# K5_2: 사업 (lower value indicates more related to men)
# K5_3: 육아 (lower value indicates more related to men)
# K5_4: 정치 
df <- df %>% mutate(gender_norm_1 = ifelse(K5_1==5| K5_1==6| K5_1==7, 1, 0),
                    gender_norm_2 = ifelse(K5_2==1| K5_2==2 | K5_2 ==3, 1, 0),
                    gender_norm_3 = ifelse(K5_3==5| K5_3==6| K5_3==7, 1, 0),
                    gender_norm_4 = ifelse(K5_4==1| K5_4==2| K5_4==3, 1, 0))

df$age.group <- NA
df[which(df$age < 30), "age.group"] <- " under 30"
df[which(df$age > 29 & df$age < 40), "age.group"] <- "30s"
df[which(df$age > 39 & df$age < 50), "age.group"] <- "40s"
df[which(df$age > 49 & df$age < 60), "age.group"] <- "50s"
df[which(df$age > 59 ), "age.group"] <- "60 and older"


df.men <- df %>% filter(BQ2=="1")
df.women <- df %>% filter(BQ2=="2")

#### Table 1, Panel (b) ####
model1 <- lm(as.numeric(gender_quota_att) ~ treated, df.men) 
model2 <- lm(as.numeric(gender_quota_att) ~ treated * young1, df.men) 
model3 <- lm(as.numeric(gender_quota_att) ~ treated, df.women) 
model4 <- lm(as.numeric(gender_quota_att) ~ treated * young1, df.women) 



#### Figure 3: Marginal Effects Plot ####
mod1.1 <- lm(as.numeric(gender_quota_att) ~ treated*young1 +  as.factor(wave) , csid)
summary(mod1.1)
estimate <- c(coef(mod1.1)[2], coef(mod1.1)[2]+coef(mod1.1)[5])
estimate
se <- c(sqrt(vcov(mod1.1)[2,2]), sqrt(vcov(mod1.1)[2,2]  + vcov(mod1.1)[5,5] + 2*vcov(mod1.1)[2,5]))
qt(0.95, df=930);qt(0.975, df=930)

ci.high.90 <- estimate  + 1.646*se
ci.low.90 <- estimate   - 1.646*se
ci.high.95 <- estimate + 1.963*se
ci.low.95 <- estimate - 1.963*se

coef.data.1 <- as.data.frame(cbind(estimate,  ci.low.90, ci.high.90, ci.low.95, ci.high.95))
coef.data.1$group <- c("Old", "Young")

mod1.2 <- lm(as.numeric(gender_quota_att) ~ treated *  young1, df.men) 
summary(mod1.2)

estimate <- c(coef(mod1.2)[2], coef(mod1.2)[2]+coef(mod1.2)[4])
se <- c(sqrt(vcov(mod1.2)[2,2]), sqrt(vcov(mod1.2)[2,2]  + vcov(mod1.2)[4,4] + 2*vcov(mod1.2)[2,4]))
ci.high.90 <- estimate  + 1.646*se
ci.low.90 <- estimate   - 1.646*se
ci.high.95 <- estimate + 1.963*se
ci.low.95 <- estimate - 1.963*se
coef.data.2 <- as.data.frame(cbind(estimate,  ci.low.90, ci.high.90, ci.low.95, ci.high.95))
coef.data.2$group <- c("Old", "Young")
coef.data <- rbind(coef.data.1, coef.data.2)
coef.data$study <- rep(c("Study 1, Male (N=930)", "Study 2, Male (N=478)"), each=2)

mod1.3 <- lm(as.numeric(gender_quota_att) ~ treated *  young1, df.women) 
estimate <- c(coef(mod1.3)[2], coef(mod1.3)[2]+coef(mod1.3)[4])
se <- c(sqrt(vcov(mod1.3)[2,2]), sqrt(vcov(mod1.3)[2,2] + vcov(mod1.3)[4,4] + 2*vcov(mod1.3)[2,4]))
ci.high.90 <- estimate  + 1.646*se
ci.low.90 <- estimate   - 1.646*se
ci.high.95 <- estimate + 1.963*se
ci.low.95 <- estimate - 1.963*se

coef.data.3 <- as.data.frame(cbind(estimate,  ci.low.90, ci.high.90, ci.low.95, ci.high.95))
coef.data.3$group <- c("Old", "Young")
coef.data.3$study <- rep("Study 2, Female (N=454)", 2)

coef.data <- rbind(coef.data, coef.data.3)
coef.data <- coef.data %>% mutate(study = fct_relevel(factor(study), "Study 1, Male (N=930)", "Study 2, Male (N=478)", "Study 2, Female (N=454)"))

ggplot(coef.data, aes(y=estimate, x=group, ymin = ci.low.95, ymax = ci.high.95)) + 
  facet_wrap(vars(study)) + 
  geom_pointrange(col="grey80") + 
  geom_pointrange(aes(ymin = ci.low.90, ymax=ci.high.90)) +
  geom_hline(yintercept = 0, lty=2, col="grey") + 
  theme_bw() + theme(panel.grid.major = element_blank()) + 
  labs(x = "Marginal Treatment Effect",
       y = "")

#### Table 2: Other Gender Equality Policies ####
model1 <- lm(as.numeric(equal_pay) ~ treated * young1, df.men) 
model2 <- lm(as.numeric(equal_pay) ~ treated * young1, df.women) 
model3 <- lm(as.numeric(corporate_quota) ~ treated * young1, df.men) 
model4 <- lm(as.numeric(corporate_quota) ~ treated * young1, df.women) 
stargazer(model1, model2, model3, model4,
          type = "latex", header=FALSE, title = "The Effect of Status Threat on Support for Gender Equality Policies",
          dep.var.labels = c("Equal Pay", "Corporate Quota"), star.cutoffs = c(0.05, 0.001, 0.0001))

#### Figure 2: Economic Insecurity by Education ####
df <- df %>% mutate(insecure = case_when(K9==1 ~ 4,
                                         K9==2 ~ 3,
                                         K9==3 ~ 2,
                                         K9==4 ~ 1),
                    college = case_when(education==1 ~ 0,
                                        education==2 ~ 0,
                                        education==3 ~ 1,
                                        education==4 ~ 1,
                                        education==5 ~ 1,
                                        education==6 ~ 1))
# For these figures, we are presenting responses of Control group only
df <- df[which(df$treated==0),]

p1 <- ggplot(df, aes(x=as.factor(college), insecure, fill= reorder(young1, desc(young1)))) + 
  stat_summary(fun = mean, geom="bar", width=.75, position="dodge") +
  stat_summary(fun.data = mean_cl_normal, geom= "errorbar",  width=0.1, alpha=1, size = 1, position=position_dodge(width = 0.75)) + 
  theme_bw() + scale_fill_manual(values = c("1" = "grey85", "0" = "grey55"),
                                 labels = c("Young (<40)", "Old (>=40)"),
                                 name = c("")) +  
  scale_x_discrete(labels= c("No College", "College")) +
  theme(axis.text = element_text(size=12), legend.text = element_text(size = 18), legend.position = "bottom",
        axis.text.y = element_text(size = 18), axis.title.y = element_text(size = 15)) + xlab("") + ylab("") + ggtitle("Economic Insecurity") + theme(text = element_text(size=12))+ coord_cartesian(ylim=c(0, 4))  
p1
t.test(as.numeric(insecure)~young1, subset(df, college==0))
t.test(as.numeric(insecure)~young1, subset(df, college==1))

